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^ . Abstract 



• In this paper, we present a regularization to ID Grad's moment system to achieve 

global hyperbolicity. The regularization is based on the observation that the charac- 
teristic polynomial of the Jacobian of the flux in Grad's moment system is independent 
of the intermediate coefficients in the Hcrmitc expansion. The method is not relied on 
the form of the collision at all, thus this regularization is applicable to the system with- 
out collision terms. Moreover, the proposed approach is proved to be the unique one if 
■ only the last moment equation is allowed to be alternated to match the condition that 

the characteristic speeds coincide with the Gauss-Hermitc interpolation points. The 
hyperbolic structure of the regularized system, including the signal speeds, Riemann 
_ invariants and the properties of the characteristic waves including the rarefaction wave, 

^ ■ contact discontinuity and shock are provided in the perfect formations. 

. Keywords: Grad's moment system; regularization; global hyperbolicity; charactcris- 

'nI" ■ tic wave 

■ 

1 Introduction 

Nowadays, the kinetic gas theory is drawing increasing attentions in the high-tech fields. 
The kinetic theory is considered as a mesoscopic description of fluids, which is based 
^ ■ on the classical Boltzmann equation from statistical physics. However, a full accurate 

mesoscopic model is still too complex for lots of problems. During a long period of time, 
people have been looking for a median model between the classical macroscopic equations 
and the Boltzmann equation. This can be tracked back to the work of Burnett [5]. As 
is well known, the Burnett equations are later proved to be linearly unstable by Bobylev 
[3]. Another way leading to linearly stable intermediate models is the moment method 
proposed by Grad [9]. Since this method was discarded by Grad himself, very few works 
contributed to this area in the last century. However, this field is becoming active in the 
recent years, since people find that some traditional difficulties in the moment equations 
can be ignored by some regularizations to these models, e.g. |12^ \T0\ \T9\ I23j. 

This paper is focusing on a major accusation against the moment method — the lack of 
global hyperbolicity for Grad's moment system. This deficiency directly causes blow-ups 
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when the distribution is far away from the equihbrium state. It has been reported that 
increasing the number of moments shows no improvements in the numerical experiments 
[8]. Levermore's work |12j gave a theoretical way to the general globally hyperbolic moment 
equations, while it is still far from practical use due to the lack of an analytical form of his 
model. Later, using the Pearson- Type-IV distribution, Torrilhon [23] also proposed a 13- 
moment system, which is globally hyperbolic when reduced to the one-dimensional case, 
but its generalization to large number moment systems seems to be difficult. In this work, 
we concentrate on the simple ID case and achieve a globally hyperbolic regularization to 
Grad's moment system successfully. 

The first essential observation is that the characteristic polynomial of Jacobian of the 
flux of a general Grad's moment system has a simple expression, which only depends on 
the macroscopic velocity, temperature, and two other coefficients in the Hermite expan- 
sion of highest orders. This amazing result directly leads to the possibility of a globally 
hyperbolic regularization. It is found that these two coefficients take the eigenvalues away 
from the real axis, resulting in the non-hyperbolicity. We discover an elegant modifica- 
tion to the last equation of the moment system to eliminate the terms involving these 
two terms in the characteristic polynomial and obtain a globally hyperbolic system. This 
new hyperbolic system has lots of fascinating properties. All characteristic fields are ei- 
ther genuinely nonlinear or linearly degenerate. The investigation into the three kinds 
of elementary characteristic waves (rarefaction waves, contact discontinuities, and shock 
waves) illustrates substantial similarities with Euler equations. The regularization pro- 
posed is very different from the classical way which tries to give a reasonable recovery of 
the truncated moments, which is justified in the view of characteristic speeds and order of 
accuracy. The convergence in the number of moments is illustrated through the numerical 
study of a shock tube problem. 

The rest of this paper is arranged as follows: in Section [21 the Boltzmann equation and 
the moment method are revised. In Section [3l a detailed investigation on the hyperbolicity 
of ID Grad's moment system is carried out. The regularization of the ID Grad's moment 
system to achieve global hyperbolicity is derived in Section [H with detailed discussion on 
its properties. A short discussion on the moment equations with collision terms is put 
forward in Section [SJ Section [S] is devoted to the numerical study of a shock tube problem. 
Finally, some concluding remarks are given in Section [71 

2 The moment method in kinetic theory 

In the kinetic gas theory, the state of a gas on the microscopic level is described by the 
velocity distribution function on each spatial point a; G C M.^. For a time-evolving 
problem, the distribution function can be described as 

F:IR+ xJlxR^^E+U{0}, {t,x,^) ^ F{t,x,$,), (2.1) 

where t is the time and ^ denotes the velocity of microscopic gas particles. As in [9], we 
introduce the mass density 

f{t,x,^) = mF{t,x,^), (2.2) 

m being the mass of the molecule. The physical case is D = 3, while in this paper, we 
only consider a ID model problem with D = 1. Thus, x and $, will be written in plain 
font as X and later on. 
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2.1 The Boltzmann equation and conservation laws 

The mass density / satisfies the Boltzmann equation, which reads 

f+«| = Q(/,/), (2,3) 

where Q{f, f) is the cohision term with a complex expression, which models the interaction 
between particles. In most part of this paper, we only consider the collisionless case, thus 
Qifj /) = will be assumed if not specified. However, the readers may keep in mind 
that our final aim is to provide an improved description of the Boltzmann equation with 
collision term using the moment method, and we will return to this topic in Section [SI 

The basic variables, including the density, the momentum density and total energy 
density, are defined as 



pit,x) = [ f{t,x,Od^, 



pit,x)u{t,x) = / ^/(t,x,e)d^, (2.4) 
lp{t,x)\u{t,x)\' + \p{t,x)9{t,x) = [ hc\^f{t,x,Od^- 

^ ^ JR 

Here u is the macroscopic velocity, and Q is the multiplication of gas constant and temper- 
ature. Multiplying the Boltzmann equation (j2.3p by (1, ,^^/2)-^, integrating both sides 
over R with respect to ^, and then making some simplifications, we get the non-conservative 
form of the conservation laws as 

'^P dp du ^ ^ \ 

,| + |E+^^=0, (2.5b) 

ot OX ox 

1 89 1 89 dq du , , 

2^5^+2^"a^ + 5^+^9^ = °' ^^-''^ 

where p is the pressure and q is the heat flux. They are deflned as 

p = p9, q=\h^-uffdt (2.6) 
^ JR 

2.2 The moment method 

The moment method was raised by Grad in [9], where a thirteen moment system was 
introduced. However, systems with large moment numbers are not investigated until 
recently (e.g. [Ml El El E] ) . Here we use the notations in [6l[8], and expand the f{t,x,(,) 
as 

fit, x,0 = Yl fk{t, x)neit,.).k ( ^~^M , (2.7) 
where ^ 
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where Hck is the k-th. Hermite polynomial, defined by 

He,{x) = (-l)^exp (:^) ^ exp (-^) . (2.9) 

Based on this expansion, some simple properties can be deduced: 

fo = P, /i = /2 = 0, g = 3/3. (2.10) 

If we put (|2.7p into the Boltzmann equation (|2.3p . the equation for each moment can be 
deduced as 

-m - ^'-'-pd^ + + '^^'d^ + U ^ J d-x ^2 11) 

p ox ox ox Ox 

For details, we refer the readers to [8]. The conservation laws ()2.5p together with ()2.1ip 
form a moment system with infinite number of equations. In order to get a closed system 
with finite number of equations, one can follow Grad's idea [9] and let Jm+i = for some 
M ^ 3. Thus a closed system with M + 1 moments is obtained. 



3 Hyperbolicity of Grad's moment systems 

A ID quasilinear system 

|.A(,)|.„ ,3.1) 

is pronounced to be hyperbolic for a particular if the matrix A(qQ) is diagonalizable 
with real eigenvalues. For Grad's systems, the hyperbolicity can only be obtained in the 
vicinity of Maxwellian [151 HI 123] ■ The loss of hyperbolicity makes the Grad's system 
overdetermined for strongly non-equilibrium gases, and severely restrict the application of 
moment methods. In this section, we are going to study the ID model problem and find 
the way in which high order moments affect the hyperbolicity of moment system. 

Let WM = {p,u,e,f3,--- JmV e M G N and M ^ 2. The Grad's moment 

system (|2.5p and (|2.1ip with fu+i = is then written as 

where Am is a lower Hessenberg matrix defined in (j3.5p . We write the matrix in a 
simplified formation with a translation and similarity transformation. Let us define 

A = diag|l,pri/2,ip^-i,0-3/2,... ,0-^V2|^ 5i = -4?' i = 3,---,M. (3.3) 
Then 



,2/-- ' /' pen 



Am = nI + \/^A-^AMA, (3.4) 
where Am is defined in (|3.6p . Thus, if 

Aj, i = l,---,M + l 
are all the eigenvalues of Am , then 

u + XjVe, j = !,■■■ ,M + 1 

are all the eigenvalues of Am • 
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u p 

e/p u 1 

26* u 6/p 

4/3 pe/2 u 4 

-ofs/p 5/4 3/3/2 e n 5 .... 

efM-2/p MfM-1 i[(M-2)/A/„2 + e/A/-4] -3/M-3/P 9 u M 

OfM-i/p (M + 1)/m i[(M-l)/M-i + e/A./-3] -3/a/_2/p 6 u 



\ 



J 



(3.5) 



1 

1 2 

1 3 

453 1 4 

■53 554 3^3 1 5 .... 



(3.6) 



-gM-2 MgM~i (M - 2)5M-2 + 5M-4 -35M-3 ••• 1 M 
^-gM-i (M + l)(7Af (M - l)(7M-i + 5A/-3 -^9ai-2 1 0^ 



The matrix Am can be considered as "simple" in a sense. It contains only dimension- 
less variables g3,--- ,gM with linear dependence. The diagonal elements of Am are all 
vanished, and the subdiagonal entries are all 1. The superdiagonal elements are equal to 
their row numbers. Meanwhile, apart from the tridiagonal part, only the first four columns 
are nonzero. These formation give us possibility to study its eigenvalues. 

We first present the main result of this section in Theorem [TJ In this paper, | • | is used 
to denote the determinant of a matrix. 



Theorem 1. The characteristic polynomial of Am is 

1 



XI- A 



M 



HcM+iW - i;{M + 1)1 ■ [{X' - l)5A/-i + 2Xg 



M 



(3.7) 



The result is incredibly simple, and therefore gives us a realistic possibility to make 
some kind of regularization to gain global hyperbolicity, which will be discussed in the 
next section. To proof this theorem, we need the following two lemmas. 

Lemma 1. Suppose that a square matrix A = {aij) depends on N variables xi, ■ ■ ■ ,xi\f. 
Then the partial derivatives of \A\ can be calculated as 



d\A\ 
dxk 



E 



OXjg 



1, 



,iV. 



(3.8) 



where A^^ is the {i,j)-th minor of matrix A, which is defined to be the determinant of the 
submatrix obtained by removing from A its i-th row and j-th column. 

This is a familiar result in linear algebra, and will not be proved here. 

Lemma 2. Define tridiagonal matrices 



fx 







-(J + 1) 

A 

-1 





-(i + 2) 

A 





■(i + 3) 



\ 














-1 A 



-M 
A 



^ i < M. (3.9) 



The following relations for the determinants of Dj hold: 



ID, 



A|D,+i| - {j + 1)|D,+2|, 



Proof. For ^ j ^ A/ — 1, Dj can be written as 



A -(j + l)ef 
-ei Dj+i 



(3.10) 



(3.11) 



where ei is the unit vector (1, 0, • • • , 0)"^. When A 7^ 0, since 



' I 





ei I 



A 

-ei 



D 



'i+i 



A -(j + l)e{ 

Dj+i - (i + l)A-ieieJ 



T 1 



(3.12) 
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the equality 



ID, 



A|Dj+i-(j + l)A-^eief 



(3.13) 



is obtained by taking determinants on both sides of ()3.12p . When ^ j ^ M — 2, we use 
p.lip again and get 



ID 



J I - A|Dj+i-(j + l)A~^eie{ 



A 



A-(j + l)A-i -(j + 2)ef 



-ei 



D 



J+2 



A 



A 



-(i + 2)e^ 



T 



-ei 



D 



'j+2 



+ 







-ei 



D 



i+2 



= A|D,-+i|-(j + l)|Dj+2|. 

(3.14) 

If A = 0, the continuity of |Dj| with respect to A gives the same result. □ 

Now we prove Theorem [TJ 

Proof of TheoremUl We start the proof by calculating d\Xl — Ai[j\/dgj for 3 ^ j ^ M — 3. 
From p.6p . one may find that gj only appears in five entries of the matrix. Their positions 
are 

(j + 2,1), (j + 1,2), (i + 2,3), (j+4,3), (j + 3,4), 



which are illustrated in Figure 1(a) Thus, according to Lemma [H only five terms appear 
in the right hand side of (|3.8p . Now we will consider them one by one. Below we denote 
AI — Aj\,/ = (cij), and use C*'-' to denote the {i,j)-th. minor of AI — Am- 



1. As in Figure 1(b) , C-'"'"^'^ is presented as the product of the determinants of two ma- 
trices. One is a lower triangular matrix whose diagonal elements are — 1, • • • , — (j+1), 
and the other is a lower right block of AI — Am, which is actually Dj+2 defined in 
Therefore, we obtain 



C^+2-i = (_i)^+i(j- + l)!.|D,+2|. (3.15) 
Since Cj+2,1 = gj, one has 

(_l),+2+l^|i2^^,+2,l ^ . 1 . + 1), . |D,.^2| = ij + 1)! • |D,+2|. 

(3.16) 



2. Figure 1(c) shows that C^~^^''^ is factorized into three parts: the first part is A, the 
second is a lower triangular matrix with diagonal elements — 2, • • • , — j, and the third 
one is Dj+i. Since Cj+1^2 = — (j + l)5ji we get 

(-i)^-+i+2%i^c^-+i'2 = (-1)^'+^ • [-(J + 1)] • {-ly-'jix ■ iD,+ii 

= -(i + l)!-A|D,+i|. 

3. C'-'+2'^ is illustrated in Figure 1(d), from which one finds C'-'+2'^ is the product of 
the determinants of three matrices. The first matrix is a 2 x 2 upper left block of 
AI — Am, for which we have 



A -1 
-1 A 



A^ 



1. 



(3.18) 
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And the other two blocks are similar as the last case. Using Cj+2,3 = —jdj, we have 

= i-iy^' ■ i-j) ■ i-iy^'lij + i)!(A^ - 1) • |D,+2i 



99j 



-(i + l)!-(A2-l)|D,-+2| 



(3.19) 



4. The structure of (7-^+^''^ is plotted in Figure 1(e), which is very similar as C-^^^'^. 
Therefore we directly write the result: 



-^y7+4+3 '9^^+4,3 ^■,+4.3 



-1) 



.+i.(_l).(_l)i+il(j+3)!(A2-l).|D,+4| 



-(i + 3)!-(A2-l)|D,+4|, 



(3.20) 



where we have used Cj+4^3 = —{j + 2)9j+2 — Qj- Note that we define |Dm+i| = 1 in 
order that ()3.20p is correct for j = M — 3. 

5. Similar as C-'+^'3 and C-'"'"^'^, the minor C^'^^''^ is also factorized into the determi- 
nants of three matrices as in Figure 1(f), while the first matrix is the 3x3 upper 

(3.21) 



left block of AI — Aj\/, whose determinant is 



A -1 
-1 A -2 
-1 A 



A^ - 3A. 



Thus the last term becomes 



^_^^,+3+4 5|±3^^,+3,4 ^ . 3 . (_1).-1 + 2)!(A3 - 3A) • |D,+3 



dg 



-(j + 2)!.(A3-3A)|D,+3| 



(3.22) 



Collecting (fXTH]) . (|3T7)l . (ITOD . (fM]) and IK22\i . we finally get 
9|AI-Am 



(j + 1) 



1 . 



D,+2|-A|D,+i|-|(A2-1)|D,h_2| 



(j+3)(j+2) ^^2 _ + J + 2(^3 _ 3^)|D^.^3| 



(3.23) 



This expression will be further simplified using Lemma [2j Since (|3.1ip also holds for 
J = M — 1 if we define |Da./+i| = 1, the following relation is deduced: 



|D,+2|-A|D,+i|-|(A2-1)|D,-+2| 



3 rx2 



|D,-+2| - A(A|D,+2| - (j + 2)|D,+3|) - ^(A' - l)|Dj+2| 



i + 2,,2 



2 (A^-l)|D,+2| + (j + 2)A|D,+3| 
^(A2 - 1)[A|D,+3| - (i + 3)|D,+4|] + (i + 2)A|D,-+3| 
^(A3 - 3A)|D,,3| + (1±M±^(A^ _ i)|D,,4|. 



(3.24) 
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Substituting this equation into ()3.23p , we conclude 



0|AI-Am| ^ 3^j^M-3. (3.25) 

It is clear that 53, • • • ,gM-3 do not appear in the characteristic polynomial of Am- 

For j = M — 2, M — 1, M, the entries containing gj still locate in the matrix as Figure 
1(a) , while some items are missing due to the cut-off. Therefore, if we define |Dj| = for 



i > M + 1, then (f3:23|) still applies for j = M - 2,M - 1,M. Note that such definition 
leads to 

|Dm| = A|Dm+i| - (M + 1)|Dm+2|, (3.26) 
therefore gM-2 does not appear in |AI — Am\ either. Moreover, we have 

d\XI - Am 



dgM~ 



Ml 



Dm+i| - A|Dm| - ^^-^{X"^ - l)|Dj+2| 



= -(^(A^-l 

d\XI-AM\ 

dgM 



(3.27) 

(M + 1)! • (-A|Dm+i|) = -(M + 1)! • A. (3.28) 



Since ()3.27p and (|3.28p hold for any gj, 3 < j < M, we write the characteristic polynomial 
of Aj\/ as 

|AI - AmI = C(A) - (^ + 1)' [(a2 _ + 2XgM], (3.29) 

where C(A) is a function of A. 

Now it only remains to determine C(A), which is done by assigning (73, • • • ,gM to be 
zero, and then calculating the characteristic polynomial of Am- In this case, it is easy to 
find 

|AI- AmI =C(A) = |Do|, if ^3 = ... = 5Af = 0. (3.30) 
Meanwhile, the following relation between |Dj| and Hermite polynomials is discovered: 

|Dm+i| = HeoiX) = 1, |Dm| = Hei{X) = A, 

|D,-| = A|Dj+i| -(j + l)|Dj+2|, (3.31) 
He,{X) = XHej.iiX) - {j - l)Hej.2{X). 

This reveals that 

|Dj| = HeM+i~j{X), ^ j ^ M + 1. (3.32) 

Hence C(A) = |Do| = HeM+i{X)- This completes the proof of Theorem[TJ □ 

Theorem [1] reveals that the hyperbolicity can only be obtained in a particular region 
{gM-i, gn) £ for the (M + l)-moment sytem. Since the roots of Hermite polynomials 
are all real, the origin must lie in f^M. The region Qm for M = 4 to 9 are plotted in Figure 
El among which the result for M = 4 has been obtained in |23j . agreeing with ours with 
proper scaling and translation. 

As a reference, the following corollary gives the characteristic polynomial of the original 
matrix Am- 
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• • • 

• • • • 



row j+l • ^ • 



(a) Entries containing gj 



row j+l 




(b) The (j + 2, l)-th minor of AI - A], 



row 7 + 1 < ^ ^ 

<5> 



(c) The (j + l, 2)-th minor of AI - Aj, 



row J + l • <^ 



(d) The (j + 2, 3)-th minor of AI - A], 



W J + l • <J> 



(e) The (j + 4, 3)-th minor of AI - Aa 



row J + l 




(f) The (j + 3, 4)-th mino: 



Figure 1: The Hessenberg matrix Am and its minors of the elements containing gj. 
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-0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 -0.08 -0.06 -0.04 -0.02 0.02 0.04 0.06 0.08 



(a) M = 4 



(b) M = 5 





(c) M = 6 



(d) M = 7 



xicr' 





-6^-2 2 4 6 



(e) M = 8 



(f) M = 9 



Figure 2: Hyperbolicity region of Grad's (M + l)-moment system. The x-axis is gu-i 
and the y-axis is qm- 
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Corollary 1. The characteristic polynomial of Am is 



. I X — u 

2 HCM+I 



(M + 1)! 



Proof. This can be shown by direct calculation: 



[((A-7x)2-0)/j,,_i + 2(A-^x)/m] 



(3.33) 



lAI- A 



, M+l 
I 2 



Ml = 

A — II 



AI - (nl + V^A^i Ajv/A) 

A — M 



I- A-^A^fA 



I 2 



I- A 



. I I \ — U 

2 <^ HCM+I 



(3.34) 



(M + l)! 



(A-«)2 /a./^i ^ 2(A-n) /m 



, M+l I \ — U 

2 i?eA/+i 



{M + l)\ 
2p 



pQ{M ~l)/2 ^ pQM/2 



[(^{^X-uf-e)fM^l + 2{\-u)fM] 



□ 



4 Hyperbolic moment system 

The loss of global hyperbolicity of Grad's moment system has long been considered as 
a failure of moment method. Recently, some encouraging progresses are made in this 
direction [12^ I23j . However, in the case that the number of moments is greater than 10, 
Levermore's method leads to great difficulties for the numerical implementation, since the 
moments cannot be analytically solved from the Lagrange multiplier^; and it has been 
demonstrated by Junk [TT] that the domain of definition for a realizable distribution is not 
convex. Torrilhon's method mainly focuses on 13-moment case in one space dimension, 
which seems not trivial to be extended to the general cases. To the best of our knowledge, 
no results for general moment system have been published. 

In this section, we provide the method to regularize moment system based on the 
results in Section [3] to achieve global hyberbolicity. We discuss only ID case here and the 
multi-dimensional problems will be reported soon in later papers. 



4.1 Construction of hyperbolic moment system 

For an (M + l)-moment system containing quantities {p,u,6, f^, - ■ ■ ,/a/}, the Grad's 
moment system gives accurate evolution equations for most variables expect for /m , since 
/m+1 appears in the accurate equation of /m, and is forced to be zero in Grad's closure. 
Almost all the regularization methods in references are focused on the reconstruction of 
trying to express /m+i as a function of the M + l known variables in some possible 
ways such as Chapman-Enskog expansion or realizing a positive distribution [121 [191 [8l [23] . 
In this paper, we also limit our regularization to the modification of the equation of /m- 
However, since /a/+i exists in this equation only in the form of its derivative, here we 
directly substitute dfM+i/dx with some other expression to gain global hyperbolicity. 



local system is required to be solved by Newton iteration on each grid for every time step. We refer 
the readers to [20] for details. There is no report indicating that such a system has a fast solver. 
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Corollary [T] shows that the characteristic polynomial of Aj\/ is independent of /a, • • • , 
fM-2, and its dependence of /m-i and /m can be regarded as the result of truncation. 
That is, if a Grad's system with M + 3 or more variables is considered, then /m-i and 
/m do not affect the characteristic polynomial, either. Thus, it is reasonable to modify 
the matrix Am such that its characteristic polynomial is a function only of u and 0. More 
precisely, the characteristic polynomial of the modified matrix should always be 



M+l 

e 2 HtM+i 



X — u 

~7f 



(4.1) 



which is obtained by substituting /m-i = fhi = into (j3.33p . Recalling that only the 
equation of Jm is allowed to be changed, we summarize all the requirements and raise the 
following problem: 

Find M + l functions aj = aj{w m) , j = 1, • " " + 1, such that 

M+l 



AI- A 



M 



M+l 

I 2 He M+l 



X — u 

~7f 



^P,u,9,f3,--- Jm, 



where Ejj denotes the matrix e^ej, and ej is the unit vector whose j-th com- 
ponent is equal to 1. 

If Oj = aj{wM), J = 1) ■ ■ ■ ) M+l is the solution of this problem, then a globally hyperbolic 
system can be obtained by substituting the matrix Am in (|3.2|) with 



M+l 

Am '■= Am + ^ ojEim+i,^ 



(4.2) 



The rest part of this section will be devoted to tackling this problem. 

In order to simplify the notation, we use S^'^ to denote the (i,j)-th minor of the 
matrix AI — Aj\,/, and define S{k) as its k-th order leading principal minor, which is the 
determinant of the upper-left part of AI — A^ with k rows and k columns. According to 
the expression of A m (|3.5p , it is not difficult to find 



gM+l,j ^ 



,M!, 



(-l)"-'^(A-.), 



1, 

Ml 



p 



rSU - 1) 



J 



(j - 1)! 

Now we expand the characteristic polynomial of the matrix ()4.2p as 



,M + 1. 



AI- A 



M 



M+l 

^ ajEM+ij 



M+l 

|AI - AmI - Yl {-l)^'^'^'ajS^'+''\ 



(4.3a) 
(4.3b) 
(4.3c) 

(4.4) 



In order that the above expression equals to (|4.ip . according to Corollary [H we may choose 
Gj such that 

(M + l)! 



M+l 



2p 



[((A - - d) fM-i + 2(A - n)fM] + {-^)^'^'"'cijS^'+''^ ^ 0. (4.5) 
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The leading principal minor S{k) is a polynomial in A of degree k, since it is the charac- 
teristic polynomial of the k x k upper-left block of Am- Hence, S^'^~^^'^ is a polynomial in 
A of degree j — 1, which can be observed from ()4.3p . Such observation directly leads to 



aj = 0, j = 4,---,M + l, (4.6) 

since the first term in (|4.5p is a quadratic polynomial in A. Then, we put (|4.3ap and (|4.3bp 
into ()4.5p and some simplification gives 



[(A -uf-9] (^^i±lfM^i + 03) +{X-u) [{M + 1)/m + 02] +ai = 0. (4.7) 
Now, the choices of ai, 02 and 03 are naturally given as 

ai = 0, a2 = -(M + l)/M, a3 = -^^^fM^i. (4.8) 
For simplicity, the notation TZm is introduced as follows: 

Definition 1. The regularization term based on the characteristic speed correction is de- 
noted as 

'^M^^f^A.I^./M-.I^V (4.0) 



I \ OX ox 
Then we have the following theorem: 
Theorem 2. The moment system 



dWM , . dWM rr, n (A ^^^\ 

nueM+i = (4.10) 



dt dx 

is strictly hyperbolic if > 0, and its characteristic speeds are 



u 



+ CjVe, j = l,---,M + l, (4.11) 



where cj is thej-th root of HeM+i{x). 

Proof. The equations (|4.10p can be rewritten as 

^+j^A„/f <.,EM,,,j^=0, (4.12) 

where Oj, j = 1, • • • ,M + 1 are defined in (j4.6p and (|4.8p . As we have discussed above, 
()4.ip gives the characteristic polynomial of the matrix in the parentheses, which will be 
denoted by Aj\/ below as in (|4.2p . If > 0, one has 

|AI - AmI = Q—Eeu^x (^^^ j = 0—HeM+i{c,) = 0. (4.13) 

Therefore, ()4.1ip gives all eigenvalues of Am- Since the Hermite polynomial i/ej\/+i(x) has 
M + 1 different zeros in M |17j . all cj's are distinct. Thus, the matrix Am has no duplicate 
eigenvalues, hence is diagonalizable. This indicates that (j4.10p is a strictly hyperbolic 
system. □ 
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Comparing with the exact moment system ()2.1ip , the hyperbohc system ()4.10p replaces 
dfM+i/dx by 

This is a totahy new way to regularize Grad's moment system. 

Remark 1. By modifying the last row of the matrix A^, the characteristic speeds can be 
appointed. Our regularization ()4.10p selects a special set of characteristic speeds ()4.1ip 
such that they coincide with the Gauss-Hermite interpolation points. As discussed in 
[22j . the characteristic speeds can be viewed as a sort of discretization of the distribution 
function. Therefore, the system (j4.10p is similar as the "shifted and scaled discrete velocity 
model" , with the expectation of spectral convergence when M goes to infinity. Meanwhile, 
unlike the ordinary discrete velocity model, the nonlinearity of Grad's moment systems 
introduced by shifting and scaling of the basis functions is preserved. Additionally, such 
regularization is only a slight modification based on the original Grad's moment system, 
and we will find in the next subsection that a number of interesting properties can be 
obtained. 

4.2 Characteristic waves of hyperbolic moment system 

In this part, we will focus on the Riemann problem of ()4.10p . First, we claim that all 
characteristic fields of (|4.10p is either genuinely nonlinear or linearly degenerate. To verify 
this, we write the right eigenvectors of Am in the following theorem: 



Theorem 3. The right eigenvector of Am with eigenvalue u + Cj^/9 i 



IS 



r 



(rj,!,--- ,rj-,M+i)^, j = l,---,M + l, (4.15) 



where cj is the j-th root of Hermite polynomial HeM+i{x), and rj^j. is defined as 

rj,i = P, rj^2 = CjVe, rj^3 = {c^j -1)9, 
Hek-i(ci) k-i c2 - 1 (4.16) 



{k-l)\ 2 ' 

Proof. To prove this theorem, we need only to prove 

AMVj = {u + CjV9)rj. (4.17) 

Split Am by row as Aj\/ = {aj , a^, - ■ ■ , a^^j^-^)'^ , where a^. is the /c-th row oi Am, k = 
1, 2, • • • , M + 1. Thus (jiTzP can be written as 

{airj,a2rj,- ■ ■ ,aM+irj)'^ = {u + CjV9)rj, forj = l,2,--- ,M + 1. (4.18) 

With the expression of Am, the first four rows of (|4.18p can be verified directly: 

airj = urj^i + prj^2 = p{u + CjVO) = rj^i{u + CjVO), 

a2rj = 0/p- rj^i + urj^2 + rj^s = CjVd{u + CjVO) = rj^2{u + CjVO), 

asrj = 29rj^2 + urj^a + 6/p ■ rj^^ = (c| - l)e{u + cjVe) = rj^^iu + cjVe), (4.19) 

a^rj = 4:f3rj^2 + pO/2- rj^s + urj^^ + 4rj,5 

= He3{cj)pe^/'^/6 ■ {u + CjVe) = rj- 4(n + ajVO), (only when M ^ 4). 
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For 5 < A; < M, 

dkrj = - ^^^Tj^i + kfk-irj^2 + 2)/fc_2 + 6'/fc-4]?^i,3 
^rj- 4 + Orjj.i + urjj + krjj+i. 

Then, we substitute ([H^ into (g^Q]), and get 

_ Hek-2{cj) k/2 , gefc„i(cj) , Hckjcj) 

_ gefc„i(cj) , ^ gefc-i(cj) (4.21) 

(A;-l)! ^^^^^ (A;-l)! 

+ {-Cje^l^){c,V9 + u)h-2 + [-0{cj - l)(n + CjVe)/2]h-3 
= {u + CjVO)rj^k- 

For k = M + 1, the situation is similar as 5 < /c < A/. We expand aM+ifj as 



_ ffeM-l(Cj) (A/+l)/2 geA./(Cj) M/2 

+ {-Cj9^/^){c,V9 + u)fM~i + [-e{c] - l){u + CjVe)/2]fM^2. 



(4.22) 



Here, Cj is the j-th root of Hermite polynomial HcM+iix). Hence, the recursion relation 
of Hermite polynomials gives 

HeM-i[Cj) - — . [4.26) 

Substituting this equation into ()4.22p . we get 

cLM+irj = (u + CjVe)r,j^M+i- (4.24) 

Collecting ()4.19p . (|4.21|) and (|4.24p . we finally arrive at (|4.18p . This completes the proof 
of the theorem. □ 

Corollary 2. iJac/i characteristic field of the hyperbolic system (|4.10p is either genuinely 
nonlinear or linearly degenerate. 

Proof. Let Sj = u + Cj\fQ., and we only need to verify that either Vw]^,Sj ■ rj = or 
^■wi{Sj ' fj ^ ^ holds. Since 

V^M^i= (o,l,^^,0,--- ,o) , (4.25) 

we have 

V^„s, • r, = c,Ve + ic,(c2 - 1) = ^c,(c| + 1)^/^. (4.26) 

If Cj is zero, the right hand side vanishes, while if Cj is nonzero, it is clear that ^ WM^j''''j ^ 
0. □ 
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This corollary indicates the simplicity of characteristic waves in the solution of Riemann 
problems. Consider the following Riemann problem: 



dWM , i dWM 

+ Mi- 



di 



dx 



w 



A/' 



0, 



X <0, 



(4.27) 



A typical solution of this problem is the composition of at most M + 2 intermediate states 



w 



.N 

M: 



W 



N+1 
M 



W 



R 
A/' 



iV ^ M, 



which are connected by elementary waves: rarefaction waves, contact discontinuities, 
or shock waves. In order to get a full understanding of the hyperbolic moment system, 
these waves will be studied respectively below. 



4.2.1 Rarefaction waves 

As all hyperbolic systems, the integral curves and the Riemann invariants are the major 
objects for the investigation of rarefaction waves. The parameterization of an integral 
curve of the vector field satisfies 



w 



(4.28) 



where C is the parameter, and 
wm{C) = 



Pio,m,m,f3io,--- j~Mic) 



T 



denotes the integral curve in the (M+l)-dimensional phase space. For a given point w^,^ in 
the phase space, the integral curve through w^^ can actually be analytically solved. Here 
we do not intend to write down the complete expressions, while the analytical solutions 
of /o(C)) u{Q and 6{C) are given as 



/5(C) = pOexp(C), 



2c.- 



exp 




ef(C) = 0°exp((c|-l)C) 



(4.29a) 
(4.29b) 
(4.29c) 



It is easy to verify that (j4.29p satisfies the first three equations of (j4.28p . Note that in 
()4.28p . only p, 6 and fj-2, fj-i appear in the right hand side of /j's equation, j = 3, • • • , M. 
Therefore, if the complete solution of wi[j{() is needed, one can solve fj{C) by explicit 
integration. Now we use ()4.29p to give the j-th eigenvalue of Aj[f {w j[f {()) as 



n C2 + 1 

cj — I 



exp 



-c 



. (4.30) 



It is not difficult to prove that Sj{'WM{C)) ^ *i(''^A/) ^^"^ only if CjC, ^ 0, which is 
helpful to predicate which part of the integral curve satisfies the entropy condition. And 
substitution of (I4.29bp into dOOll gives 



SjiwAliO) - Sj{w 



M) 



(4.31) 
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Hence, Sj{wM{C)) ^ ^ji'^^i) holds if and only if u{() ^ u^. Therefore, if the left state 
wj^ and the right state w^^ are connected by a single rarefaction wave, < has to 
be satisfied, since the entropy condition requires Sj{w^,^) < Sj{w^). Now let us turn to 
the pressure p. Equations (j4.29ap and (j4.29cp show that 

p(C) = p(C)e~(C)=/exp(c2C). (4.32) 

Therefore, the pressures on both sides of a rarefaction wave should satisfy 

pL ^ ^£ ^. y Q. 

> p^, if Cj < 0. 
Here we point out that the sign of Cj is as 

>0, ifi>(M + l)/2, 

= 0, ifi = (M + l)/2, (4.33) 
<0, ifi<(M + l)/2. 

It is interesting that Riemann invariants exist for all genuinely nonlinear fields, and 
the following theorem gives its expressions. 

Theorem 4. For hyperbolic moment system (|4.1U|) . the Riemann invariants for the j- 
family are 



k 

Rk = Ckfipe''/^ + J2 Ck4^^^''-'^^^ k = 3,--- ,M. 



(4.34) 



1=3 

where i is defined recursively as 



2cj 



Ck,k — 1, Ck,k-i — -2 7! (4.35a) 



cj-l 



Ck,i = ^ \ CkA^2 + Cfc,i+i j , i = 3, • • • , A; - 2, (4.35b) 

C'.o=^rT|^E^C.,,.. ,4,35c) 

^ 3' i=3 

Proof. We only need to prove 

V^,,Rk-rj = 0, yk = l,---,M. (4.36) 
The verification in the cases k = 1 and /c = 2 is straightforward: 

V^,,i?i • r, = ^-V(S^-i) . p _ 1 p^-i/Cs^-D-i . _ 1)9 = 0, (4.37) 

V^,,ii2 ■r, = l.c,V9- ^ • {cj - 1)6 = 0. (4.38) 

(c| - 1)V0 
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If /c ^ 3, the gradient of Rk is 



fc-i 



i=3 



(4.39) 



With some rearrangement, Vw^Rk ■ fj is simplified as 

k 



Hei{cj) 



k-2 
1=3 



c?-l 



d-1 



— (A; - i)Cfc^j — CjCk,i+i ^-r — Ck,i+2 



'c?-l 



+ 



-Ck,k-1 — CjCk,k fk-l- 



f.g{k~i)/2 (4 40) 



We have that 

• (|4.35ap indicates that the last line of (|4.40p is zero; 

• ()4.35bp indicates that the second line of (|4.4U|) is zero; 

• ()4.35cp indicates that the first line of ()4.40p is zero. 
Thus (|4.36p is proved. 

4.2.2 Contact discontinuities 



□ 



According to the proof of Corollary [21 the contact discontinuities can only be found in the 
case of Cj = 0. Thus, if M is odd, no contact discontinuities exist in the characteristic 
waves. For contact discontinuities, the discussion on integral curves and Riemann invari- 
ants above is still valid. If we substitute Cj = into (|4.34p . u, p and /s can be found to 
be invariant acrossing the contact discontinuity. 

4.2.3 Shock v^aves 

Discussion of the shock waves requires additional scrupulosity. As been well known, the 
jump condition on the shock wave is sensitive to the form of the hyperbolic equations. 
Therefore, before we give the Rankine-Hugoniot condition, it is necessary to rewrite ()4.10p 
in an appropriate form. Though a conservative form is desired, the whole system can no 
longer be written as a conservation law since two terms are added to the last equation. 
Nevertheless, the conservative form of the first M equations remains. Thus ()4.10p can 
actually be reformulated by M conservation laws and a single non-conservative equation. 
Precisely, if we let 



{qo, ■ ■ ■ ,qM)'^, 



[ ev(e)de, j = o,-- - ,M, 

j' Jr 



(4.41) 
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(|4.10p is reformulated as 



dt dx 
The relation between q and wm is 



(4.42) 



/, = n = gi/go, ^ = 2^2/90 - te/cZo)^ (4.43) 



and -F(q) is defined as 



F(,) = (M + 1) ^'-''^ 
For convenience, we write ()4.42p in the following form: 

|.B(„|i.„. ,4.45) 

where B(q) is an (M + 1) x (M + 1) matrix. 

Since ()4.45p is still a non-conservative system, the DLM theory [13] is introduced when 
discussing the shock wave. A shock wave is a single jump discontinuity connecting two 
constant states and q^ in a genuinely nonlinear field j, and q^, q^ and the propagation 
speed of the shock wave Sj should satisfy the following conditions: 

• Generalized Rankine-Hugoniot condition: 

[5,I-B(*(i.;q^,q«))] — (z.;q^,g^)dz. = 0, (4.46) 

where I is the identity matrix of order M + 1, and ^{v] q^ , q^) is a locally Lipschitz 
mapping satisfying 

*(0;q^,g^) = q^, *(1; q^, q^) = g^. (4.47) 

We refer the readers to [13] for details. In Section [5l we will point out that the setup 
of $ is not crucial if the collision term presents. 

• Entropy condition 

5,(q^) > S, > s,(q^). (4.48) 

It is obvious that the first M rows of ()4.46p are independent of they are the same as 
the classical Rankine-Hugoniot conditions. This allows us to analyze the properties of the 
shock waves without regarding the form of 

The first and second equations of ()4.46p can be written as 

p^u^ - = Sj{p^ - p^), (4.49) 
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Since 7^ and p^u^ 7^ p^u^ (otherwise = q^), one has 

S, = - (4.51a) 

^ p\u^ f ^ pH^ - p\u^f - pH^ _ 

pLyL _ pRyR 

Putting (|4.51a|) into (|4.48|) . and multiplying both sides with (p-^ — p^)^, we get 

p^iu^ - u^){p^ - p^) > Cj{p^ - p^fVdR, (4.52a) 

p^{u^ - u^){p^ - p^) < Cj{p^ - p^fVd^. (4.52b) 

If Cj > 0, (|i32i)l gives 

{u^ - u^){p^ - p^) > 0. (4.53) 

Thus, we can divide both sides of (|4.52p by {u^ — u^){p^ — p^) without changing the 
inequality sign, and the result is: 



from which one directly has 

{pLfe^ - {p^fe^ > 0. (4.55) 

Similarly, if Cj < 0, we have 

(n^-'u^)(p^-/) <0, and {p^fe^-{p^fe^ <^. (4.56) 

Lemma 3. For hyperbolic moment system ()4.45p . if q^ and q^ are connected by a j -shock 
wave, then the following inequalities hold: 



u'^ > and { \ R , n (4-57) 

{ P < P , if Cj < 0. ^ ^ 

Proof. With some rearrangement, ()4.5ip can be reformulated as 

(/ _ /)(p^0^ - p«0«) = /p«(n^ - 7x«)2. (4.58) 

Since the right hand side of ()4.58p is positive, one but only one of the following two 
statements is true: 

1. p^ > p^ and p^O^ > p^e^; 

2. p^ < p^ and p^e^ < p^O^. 

If Cj > 0, equation ()4.55p indicates that the first statement is true. Then, we can use 
()4.53p to conclude > u^. The conclusion for the case Cj < can be proved in the same 
way. □ 

Now, we summarize all our discussions on the entropy conditions of three types of 
waves in the following theorem: 
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Theorem 5. For hyperbolic moment system ()4.10p . if the wave of the j-th family is 
elementary, then its type can he determined by the value of Cj and the macroscopic velocities 
or pressures on both sides of the wave: 





Velocity 


Pressure 


Contact discontinuity 


Cj = 0, = 


Cj = 0, p^ = p^ 


Rarefaction wave 


Cj ^0, < 


Cj > 0, p^ < p^ 


Cj <0, p^ > p^ 


Shock wave 


Cj ^0, > 


Cj > 0, p^ > p^ 


Cj <0, p^ < p^ 



Remark 2. It is not difficult to find that Euler equations are a special case of the proposed 
hyperbolic moment equations. In the case of M = 2, we have /i = /2 = thus the 
regularization vanishes. In other words, just like Grad's moment sytem, the hyperbolic 
system can be viewed as an extension of Euler equations. Actually, all the discussions 
in this section, including the eigenvalues and eigenvectors, Riemann invariants, and the 
entropy condition, are valid for the ID Euler equations with adiabatic index 7 = 3, while 
Grad's moment system is not able to preserve these criterions. In this respect, comparing 
with Grad's moment system, this regularized moment system is likely to be a more natural 
extension of Euler equations. 



5 The case with colUsion terms 

In this section, we will give a short discussion on the moment system with collision terms. 
For simplicity, the BGK collision operator [2j is considered. In this case, the Boltzmann 
equation (|2.3p becomes 

= >-/). M 

where r is the relaxation time, and /m is the Maxwellian: 

This equation leads to a very simple form of the collision term in the moment system as 

dWM , . dwM . „ 1„ 

+ Am—^ duK-MeM+i = --t'wM-, (5.3) 

where P is a diagonal matrix 

P = diag{0,0,0,l,--- ,1}, (5.4) 

and 5h = ^ corresponds to Grad's moment system, while 5// = 1 corresponds to the 
regularized moment system. Note that when considering the weak solution of ()5.3p . one 
still needs to rewrite (|5.3p as equations of q: 



^ + {j + l)^ = --V,{qo,---,q,), J=0,...,M-1, 

dq^rdFiq) 1 ^"-'^ 

"dT dx ^ M(qo, ■ ■ ■ , qki), 
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where Vj, j = 0, • • • ,M are the corresponding production terms. Then the first order 
derivative part of the last equation will still be treated using the DLM theory. 

An important index that exhibits the quality of a collisional moment system is its 
order of accuracy in term of r. The conception of "order of accuracy" is based on the 
assumption that r is a small quantity, and its precise definition can be found in |18t 17]. In 
[8], the order of magnitude for each moment has been deduced as 

/fc~0(rr^'/3l)^ (5.6) 

for the infinite moment system, which is obtained by the technique of Maxwellian iteration. 
It is easy to find that ()5.6p remains correct for the regularized moment system (equation 
(j5.3p with 5h = 1), since the order of TZm never exceeds the leading order term of /a//t"- 
However, when M = 3m + 1, m ^ 1, the order of accuracy of the moment system is 
actually reduced by 2 with presentation of the regularization terms. This fact is not 
difficult to obtain and will be reported elsewhere. In general, the order of accuracy still 
goes to infinity as M increases. 

Another issue is the choice of the path function which was introduced in (j4.46p . 
Let us restrict our discussion of its role in solving a Riemann problem of ()5.3p . First, we 
need to get some knowledge about the general behavior of the solution, referring to the 
careful study of the Riemann problem of 13-moment system in [21]. Roughly speaking, 
the solution shows a number of waves initially, then these waves are damping gradually, 
and eventually the solution tends to a smooth curve which is similar as the solution of 
Euler equations. The initial waves have no physical meanings due to the strong non- 
equilibrium which cannot be described by the moment system, while the solution gets 
close to the Boltzmann equation's solution only when the waves are fully dissipated. Later, 
this behavior is verified numerically for large number moment equations in [1], where the 
authors show that the speed of dissipation increases when the number of moments gets 
larger. It is expected that this also describes the evolution of regularized moment system. 
Based on [211 [T]. we have the following assertions for the regularized moment system: 

1. If subshocks appear in the solution, the choice of ^ indeed makes sense. In this 
situation, the system is inadequate for the description of the physical process, saying 
M needs to be increased. 

2. $ affects the solution when the time t is very small. However, such solution has no 
physical significance, either. Only when the solution gets close enough to a smooth 
function, the moment system starts to show its ability to describe physics. Note that 
the smooth solution is independent of therefore, $ only affects the way in which 
the waves are damped, but does not affect the intrinsic constituent of the solution. 

These two assertions indicate that the choice of $ is not crucial in solving a Riemann 
problem. We can simply use a linear function to connect any two states such that the 
numerical schemes can be constructed easily. 
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6 Numerical experimentation for a shock tube problem 



In this section, a shock tube problem is studied numericahy to show the behavior of the 
hyperbohc moment systems. We consider the fohowing Riemann problem: 

dWM , I dwM 1„ 

Ot ox T 

( wjij, x<0, (6-1) 
wm{0,x) = < 

[ wfi, x>0, 

where P is defined in ()5.4p and the initial left and right states are 

toi, = (7,0,l,0,--- ,0f , wf, = {1,0,1,0,- ■■ ,0f. (6.2) 

The relaxation time is chosen as r = Kn/p. Here two different cases Kn = 0.05 and 
Kn = 0.5 are considered. A nonconservative version of the HLL scheme [16] is employed 
to discretize the moment system. 

The numerical results for Kn = 0.05 with M ranging from 2 to 10 are listed in FigureEl 
in which the thin black lines are the numerical results of the hyperbolic moment equations 
(HME), and the thick gray lines are the results of Mieussens' discrete velocity model 
(DVM) [13], provided as reference solutions. The profiles of p, u and p are drawn. It 
is clear that the solutions of hyperbolic moment systems converge to the solution of the 
Boltzmann equation when M increases. Note that when M = 2, the hyperbolic moment 
system is equivalent to the Euler equations, and the contact discontinuities and the shocks 
are obvious. When M = 3, a shock can still be found near x = 0.75. When M is greater 
than 5, the discontinuities are fully damped. This agrees with Torillhon's theory [1] that 
the discontinuities are damped faster when M is larger. 

For a larger Knudsen number Kn = 0.5, the results are shown in Figure [H These 
results can also be considered as the solutions at t = 0.03 in the case of Kn = 0.05 (with 
proper scaling in the x direction). Thus these actually show the start-up phases of a shock 
tube by moment approximation. The discontinuities are clear for all choices of AI, and 
the convergence can also be observed by eyes. 



7 Concluding remarks 

We regularize the ID Grad's moment system to achieve global hyperbolicity for arbitary 
order expansion. Fully investigations to the characteristic waves show that this set of 
equations may be a natural extension of Euler equations. Actually, the approach in this 
paper has been extended to two or three dimensional Grad's moment system, and the 
result is reported in a following paper. 
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Figure 3: Numerical results of the shock tube problem for Kn = 0.05. The left y-axis is 
for p and p, and the right y-axis is for u. 
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